# ============================================================
# fsQCA SENSITIVITY ANALYSIS
# Medicaid Unwinding Procedural Disenrollment Study
# March 2026
# ============================================================
#
# BASELINE (from fsQCA_Analysis_Mar2026.R, Model A, N=50):
#   C1: crisp — 0→0.0, 1→0.499, 2→1.0
#   C2: calibrate(c2_sum,   3,    4.5,  7)
#   C3: calibrate(c3_nobb,  0.23, 0.281, 0.34)
#   C4: calibrate(c4_abandon, 0.05, 0.15, 0.35)
#   OUT: calibrate(index2_mm, 0.20, 0.50, 0.80)
#   incl.cut = 0.80 | n.cut = 1
#
# TESTS IMPLEMENTED:
#   S1. Consistency threshold: 0.75 / 0.80* / 0.85
#   S2. Frequency threshold:   n=1* / n=2
#   S3. Outcome anchors:       tight (0.15/0.50/0.85)
#                              baseline (0.20/0.50/0.80)*
#                              loose  (0.25/0.50/0.75)
#   S4. C4 crossover:          0.12 / 0.15* / 0.18
#   S5. C3 crossover:          0.27 / 0.281* / 0.30
#
# * = baseline specification
# ============================================================

library(QCA)
library(SetMethods)
library(readxl)
library(writexl)
library(dplyr)

setwd("~/Library/CloudStorage/Box-Box/_Medicaid/_Soomi_2026")
output_dir <- "Sensitivity/"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)

# ============================================================
# 1. LOAD DATA — Model A (N=50, SD excluded)
# ============================================================

df_raw <- read_excel("data_mar2026_v3.xlsx", sheet = "Data")
df <- df_raw %>% filter(id != "SD")
cat("N =", nrow(df), "\n\n")


# ============================================================
# 2. BASELINE CALIBRATION (identical to main analysis)
# ============================================================

calibrate_baseline <- function(df) {
  # C1: crisp three-value (no change across any sensitivity test)
  df$fs_c1 <- ifelse(df$c1_incomest == 0, 0.0,
                     ifelse(df$c1_incomest == 1, 0.499, 1.0))
  
  # C2: composite documentation burden (no sensitivity test — crossover
  #     is at 4.5, a value no state occupies; shifting it is a non-test)
  df$fs_c2 <- calibrate(df$c2_sum, thresholds = c(3, 4.5, 7))
  
  # C3: digital access barrier — baseline crossover 0.281
  df$fs_c3 <- calibrate(df$c3_nobb, thresholds = c(0.23, 0.281, 0.34))
  
  # C4: interaction burden — baseline crossover 0.15
  df$fs_c4 <- calibrate(df$c4_abandon, thresholds = c(0.05, 0.15, 0.35))
  
  # Outcome: min-max normalize then calibrate — baseline anchors 0.20/0.50/0.80
  df$index2_mm <- (df$index2 - min(df$index2)) / (max(df$index2) - min(df$index2))
  df$fs_out    <- calibrate(df$index2_mm, thresholds = c(0.20, 0.50, 0.80))
  
  return(df)
}

df <- calibrate_baseline(df)
rownames(df) <- df$id


# ============================================================
# 3. HELPER FUNCTIONS
# ============================================================

# Run truth table + minimization, return compact summary
run_one <- function(df, incl.cut = 0.80, n.cut = 1, label = "") {
  
  cat("\n", strrep("=", 60), "\n")
  cat(" ", label, "\n")
  cat(" incl.cut =", incl.cut, "| n.cut =", n.cut, "\n")
  cat(strrep("=", 60), "\n")
  
  # --- Positive outcome ---
  TT_pos <- tryCatch(
    truthTable(data = df, outcome = "fs_out",
               conditions = c("fs_c1","fs_c2","fs_c3","fs_c4"),
               incl.cut = incl.cut, n.cut = n.cut,
               show.cases = TRUE, complete = FALSE),
    error = function(e) { cat("  ERROR (positive TT):", conditionMessage(e), "\n"); NULL }
  )
  
  sol_pos_expr <- NA_character_
  incl_pos     <- NA_real_
  cov_pos      <- NA_real_
  pri_pos      <- NA_real_
  
  if (!is.null(TT_pos)) {
    cat("\n-- Truth Table (Positive) --\n")
    print(TT_pos)
    sol_pos <- tryCatch(minimize(TT_pos, details = TRUE, include = "?"),
                        error = function(e) { cat("  Minimization error:", conditionMessage(e), "\n"); NULL })
    if (!is.null(sol_pos)) {
      cat("\n-- Solution (Positive) --\n")
      print(sol_pos)
      # Extract solution formula as string
      sol_pos_expr <- tryCatch(
        paste(sol_pos$solution[[1]], collapse = " + "),
        error = function(e) "parse error"
      )
      # Safe fit extractor: handles single vector OR data frame (multiple models)
      safe_fit <- function(sol) {
        tryCatch({
          pof <- sol$pof
          if (!is.null(pof) && is.data.frame(pof)) {
            last <- pof[nrow(pof), ]
            return(c(inclS = as.numeric(last[["inclS"]]),
                     covS  = as.numeric(last[["covS"]]),
                     PRI   = as.numeric(last[["PRI"]])))
          }
          return(NULL)
        }, error = function(e) NULL)
      }
      
      fit <- safe_fit(sol_pos)
      if (!is.null(fit)) {
        incl_pos <- round(fit[["inclS"]], 3)
        cov_pos  <- round(fit[["covS"]],  3)
        pri_pos  <- round(fit[["PRI"]],   3)
      }
    }
  }
  
  # --- Negative outcome ---
  TT_neg <- tryCatch(
    truthTable(data = df, outcome = "fs_out",
               conditions = c("fs_c1","fs_c2","fs_c3","fs_c4"),
               incl.cut = incl.cut, n.cut = n.cut,
               show.cases = TRUE, complete = FALSE, neg.out = TRUE),
    error = function(e) { cat("  ERROR (negative TT):", conditionMessage(e), "\n"); NULL }
  )
  
  sol_neg_expr <- NA_character_
  incl_neg     <- NA_real_
  cov_neg      <- NA_real_
  pri_neg      <- NA_real_
  
  if (!is.null(TT_neg)) {
    cat("\n-- Truth Table (Negative) --\n")
    print(TT_neg)
    sol_neg <- tryCatch(minimize(TT_neg, details = TRUE, include = "?"),
                        error = function(e) { cat("  Minimization error:", conditionMessage(e), "\n"); NULL })
    if (!is.null(sol_neg)) {
      cat("\n-- Solution (Negative) --\n")
      print(sol_neg)
      flush.console()
      sol_neg_expr <- tryCatch(
        paste(sol_neg$solution[[1]], collapse = " + "),
        error = function(e) "parse error"
      )
      fit_neg <- safe_fit(sol_neg)
      if (!is.null(fit_neg)) {
        incl_neg <- round(fit_neg[["inclS"]], 3)
        cov_neg  <- round(fit_neg[["covS"]],  3)
        pri_neg  <- round(fit_neg[["PRI"]],   3)
      }
    }
  }
  
  return(data.frame(
    label        = label,
    incl_cut     = incl.cut,
    n_cut        = n.cut,
    sol_pos      = sol_pos_expr,
    inclS_pos    = incl_pos,
    covS_pos     = cov_pos,
    PRI_pos      = pri_pos,
    sol_neg      = sol_neg_expr,
    inclS_neg    = incl_neg,
    covS_neg     = cov_neg,
    PRI_neg      = pri_neg,
    stringsAsFactors = FALSE
  ))
}


# ============================================================
# 4. S1 — CONSISTENCY THRESHOLD VARIATION
#    Baseline: incl.cut = 0.80
#    Tests:    incl.cut = 0.75  (more permissive)
#              incl.cut = 0.85  (more restrictive — hardest test)
# ============================================================

cat("\n\n")
cat("############################################################\n")
cat("  S1: CONSISTENCY THRESHOLD VARIATION\n")
cat("############################################################\n")

s1_base <- run_one(df, incl.cut = 0.80, n.cut = 1, label = "S1-Baseline  | incl=0.80, n=1")
s1_low  <- run_one(df, incl.cut = 0.75, n.cut = 1, label = "S1-Low       | incl=0.75, n=1")
s1_high <- run_one(df, incl.cut = 0.85, n.cut = 1, label = "S1-High      | incl=0.85, n=1")


# ============================================================
# 5. S2 — FREQUENCY THRESHOLD VARIATION
#    Baseline: n.cut = 1
#    Test:     n.cut = 2  (drops singleton configurations)
# ============================================================

cat("\n\n")
cat("############################################################\n")
cat("  S2: FREQUENCY THRESHOLD VARIATION\n")
cat("############################################################\n")

s2_n2 <- run_one(df, incl.cut = 0.80, n.cut = 2, label = "S2-Freq=2    | incl=0.80, n=2")


# ============================================================
# 6. S3 — OUTCOME CALIBRATION ANCHOR VARIATION
#    Baseline:  calibrate(index2_mm, 0.20, 0.50, 0.80)
#    Tight:     calibrate(index2_mm, 0.15, 0.50, 0.85)  — stricter
#    Loose:     calibrate(index2_mm, 0.25, 0.50, 0.75)  — more inclusive
# ============================================================

cat("\n\n")
cat("############################################################\n")
cat("  S3: OUTCOME CALIBRATION ANCHOR VARIATION\n")
cat("############################################################\n")

# Tight anchors
df_s3_tight <- df
df_s3_tight$fs_out <- calibrate(df_s3_tight$index2_mm, thresholds = c(0.15, 0.50, 0.85))
# Check crossover
if (any(df_s3_tight$fs_out == 0.5)) {
  cat("  WARNING: crossover conflict in S3-Tight for:",
      df_s3_tight$id[df_s3_tight$fs_out == 0.5], "\n")
}
cat("S3-Tight outcome distribution:\n")
cat("  Above 0.5 (positive set):", sum(df_s3_tight$fs_out > 0.5), "\n")
cat("  Below 0.5 (negative set):", sum(df_s3_tight$fs_out < 0.5), "\n")
rownames(df_s3_tight) <- df_s3_tight$id
s3_tight <- run_one(df_s3_tight, incl.cut = 0.80, n.cut = 1,
                    label = "S3-Tight     | OUT anchors: 0.15/0.50/0.85")

# Loose anchors
df_s3_loose <- df
df_s3_loose$fs_out <- calibrate(df_s3_loose$index2_mm, thresholds = c(0.25, 0.50, 0.75))
if (any(df_s3_loose$fs_out == 0.5)) {
  cat("  WARNING: crossover conflict in S3-Loose for:",
      df_s3_loose$id[df_s3_loose$fs_out == 0.5], "\n")
}
cat("S3-Loose outcome distribution:\n")
cat("  Above 0.5 (positive set):", sum(df_s3_loose$fs_out > 0.5), "\n")
cat("  Below 0.5 (negative set):", sum(df_s3_loose$fs_out < 0.5), "\n")
rownames(df_s3_loose) <- df_s3_loose$id
s3_loose <- run_one(df_s3_loose, incl.cut = 0.80, n.cut = 1,
                    label = "S3-Loose     | OUT anchors: 0.25/0.50/0.75")


# ============================================================
# 7. S4 — C4 CROSSOVER VARIATION (interaction burden)
#    Baseline:  calibrate(c4_abandon, 0.05, 0.15, 0.35)
#    Lower XO:  calibrate(c4_abandon, 0.05, 0.12, 0.35)
#    Higher XO: calibrate(c4_abandon, 0.05, 0.18, 0.35)
# ============================================================

cat("\n\n")
cat("############################################################\n")
cat("  S4: C4 CROSSOVER VARIATION (interaction burden)\n")
cat("############################################################\n")

# Lower crossover: 0.12 — more states classified as high burden
df_s4_low <- df
df_s4_low$fs_c4 <- calibrate(df_s4_low$c4_abandon, thresholds = c(0.05, 0.12, 0.35))
if (any(df_s4_low$fs_c4 == 0.5)) {
  cat("  WARNING: crossover conflict in S4-Low C4 for:",
      df_s4_low$id[df_s4_low$fs_c4 == 0.5], "\n")
}
rownames(df_s4_low) <- df_s4_low$id
s4_low <- run_one(df_s4_low, incl.cut = 0.80, n.cut = 1,
                  label = "S4-C4xo=0.12 | C4 crossover lowered to 0.12")

# Higher crossover: 0.18 — fewer states classified as high burden
df_s4_high <- df
df_s4_high$fs_c4 <- calibrate(df_s4_high$c4_abandon, thresholds = c(0.05, 0.18, 0.35))
if (any(df_s4_high$fs_c4 == 0.5)) {
  cat("  WARNING: crossover conflict in S4-High C4 for:",
      df_s4_high$id[df_s4_high$fs_c4 == 0.5], "\n")
}
rownames(df_s4_high) <- df_s4_high$id
s4_high <- run_one(df_s4_high, incl.cut = 0.80, n.cut = 1,
                   label = "S4-C4xo=0.18 | C4 crossover raised to 0.18")


# ============================================================
# 8. S5 — C3 CROSSOVER VARIATION (digital access barrier)
#    Baseline:  calibrate(c3_nobb, 0.23, 0.281, 0.34)
#    Lower XO:  calibrate(c3_nobb, 0.23, 0.27,  0.34)
#    Higher XO: calibrate(c3_nobb, 0.23, 0.30,  0.34)
# ============================================================

cat("\n\n")
cat("############################################################\n")
cat("  S5: C3 CROSSOVER VARIATION (digital access barrier)\n")
cat("############################################################\n")

df_s5_low <- df
df_s5_low$fs_c3 <- calibrate(df_s5_low$c3_nobb, thresholds = c(0.23, 0.27, 0.34))
if (any(df_s5_low$fs_c3 == 0.5)) {
  cat("  WARNING: crossover conflict in S5-Low C3 for:",
      df_s5_low$id[df_s5_low$fs_c3 == 0.5], "\n")
}
rownames(df_s5_low) <- df_s5_low$id
s5_low <- run_one(df_s5_low, incl.cut = 0.80, n.cut = 1,
                  label = "S5-C3xo=0.27 | C3 crossover lowered to 0.27")

df_s5_high <- df
df_s5_high$fs_c3 <- calibrate(df_s5_high$c3_nobb, thresholds = c(0.23, 0.30, 0.34))
if (any(df_s5_high$fs_c3 == 0.5)) {
  cat("  WARNING: crossover conflict in S5-High C3 for:",
      df_s5_high$id[df_s5_high$fs_c3 == 0.5], "\n")
}
rownames(df_s5_high) <- df_s5_high$id
s5_high <- run_one(df_s5_high, incl.cut = 0.80, n.cut = 1,
                   label = "S5-C3xo=0.30 | C3 crossover raised to 0.30")


# ============================================================
# 9. COMBINED HARD TEST
#    Most demanding single specification:
#      incl.cut = 0.85  (strictest consistency threshold)
#      n.cut    = 2     (no singletons)
#      OUT      = tight anchors (0.15/0.50/0.85)
#      C4 xo    = 0.12  (more states classified as high burden)
#    If the robust core survives this, it is very strong.
# ============================================================

cat("\n\n")
cat("############################################################\n")
cat("  S6: COMBINED HARD TEST\n")
cat("  incl=0.85 | n=2 | OUT tight | C4 xo=0.12\n")
cat("############################################################\n")

df_hard <- df
df_hard$fs_out <- calibrate(df_hard$index2_mm, thresholds = c(0.15, 0.50, 0.85))
df_hard$fs_c4  <- calibrate(df_hard$c4_abandon, thresholds = c(0.05, 0.12, 0.35))
if (any(df_hard$fs_out == 0.5)) {
  cat("  WARNING: crossover conflict in hard test outcome for:",
      df_hard$id[df_hard$fs_out == 0.5], "\n")
}
if (any(df_hard$fs_c4 == 0.5)) {
  cat("  WARNING: crossover conflict in hard test C4 for:",
      df_hard$id[df_hard$fs_c4 == 0.5], "\n")
}
rownames(df_hard) <- df_hard$id
s6_hard <- run_one(df_hard, incl.cut = 0.85, n.cut = 2,
                   label = "S6-Hard      | incl=0.85, n=2, tight OUT, C4xo=0.12")


# ============================================================
# 10. SUMMARY TABLE
# ============================================================

cat("\n\n")
cat("############################################################\n")
cat("  SENSITIVITY ANALYSIS SUMMARY TABLE\n")
cat("############################################################\n")

summary_tbl <- bind_rows(
  s1_base,
  s1_low,
  s1_high,
  s2_n2,
  s3_tight,
  s3_loose,
  s4_low,
  s4_high,
  s5_low,
  s5_high,
  s6_hard
)

cat("\n--- Positive Outcome Solutions ---\n")
print(summary_tbl %>% select(label, incl_cut, n_cut, sol_pos, inclS_pos, covS_pos, PRI_pos),
      right = FALSE)

cat("\n--- Negative Outcome Solutions ---\n")
print(summary_tbl %>% select(label, incl_cut, n_cut, sol_neg, inclS_neg, covS_neg, PRI_neg),
      right = FALSE)

# Export summary
write_xlsx(summary_tbl, path = paste0(output_dir, "sensitivity_summary.xlsx"))
cat("\nSummary table saved to:", paste0(output_dir, "sensitivity_summary.xlsx"), "\n")


# ============================================================
# 11. CROSSOVER CONFLICT CHECK FOR ALL SPECS
# ============================================================

cat("\n\n--- Crossover Conflict Summary ---\n")
specs <- list(
  list(df=df,         label="Baseline"),
  list(df=df_s3_tight,label="S3-Tight (OUT 0.15/0.85)"),
  list(df=df_s3_loose,label="S3-Loose (OUT 0.25/0.75)"),
  list(df=df_s4_low,  label="S4-C4xo=0.12"),
  list(df=df_s4_high, label="S4-C4xo=0.18"),
  list(df=df_s5_low,  label="S5-C3xo=0.27"),
  list(df=df_s5_high, label="S5-C3xo=0.30"),
  list(df=df_hard,    label="S6-Hard")
)
vars  <- c("fs_c1","fs_c2","fs_c3","fs_c4","fs_out")
cnames <- c("C1","C2","C3","C4","OUT")

for (sp in specs) {
  conflicts <- c()
  for (i in seq_along(vars)) {
    v <- vars[i]
    if (v %in% names(sp$df)) {
      hit <- sp$df$id[sp$df[[v]] == 0.5]
      if (length(hit) > 0) conflicts <- c(conflicts, paste0(cnames[i], ": ", paste(hit, collapse=",")))
    }
  }
  if (length(conflicts) == 0) {
    cat("  ✓", sp$label, "— no conflicts\n")
  } else {
    cat("  ⚠", sp$label, "— CONFLICTS:", paste(conflicts, collapse=" | "), "\n")
  }
}

cat("\n\nSensitivity analysis complete.\n")
cat("All outputs saved to:", output_dir, "\n")

# ============================================================
# QUICKFIX: Run S3-Loose negative solution standalone
# (Run this in your R session if the main script errored before printing it)
# ============================================================
# Assumes df_s3_loose is already in your environment from S3 block above

TT_neg_loose <- truthTable(
  data       = df_s3_loose,
  outcome    = 'fs_out',
  conditions = c('fs_c1','fs_c2','fs_c3','fs_c4'),
  incl.cut   = 0.80, n.cut = 1,
  show.cases = TRUE, complete = FALSE, neg.out = TRUE
)
cat('
-- S3-Loose: Truth Table (Negative) --
')
print(TT_neg_loose)

sol_neg_loose <- minimize(TT_neg_loose, details = TRUE, include = '?')
cat('
-- S3-Loose: Solution (Negative) --
')
print(sol_neg_loose)